function [TEC, drho, tau] = NeQuickG(h1, lat1, lon1, h2, lat2, lon2, data, f)
%
% NeQuick-G Ionospheric Model
%
%DESCRIPTION:
%This function implements the computation of the TEC with NeQuickG.
%
%PROTOTYPE:
% [TEC] = NeQuickG(h1, lat1, lon1, h2, lat2, lon2, data)
% [TEC, drho, tau] = NeQuickG(h1, lat1, lon1, h2, lat2, lon2, data, f)
%
%--------------------------------------------------------------------------
% INPUTS:
%   h1         [1x1]       Point 1 Altitude          [km]  (User)
%   lat1       [1x1]       Point 1 Latitude          [deg] (User) [-90,+90]deg
%   lon1       [1x1]       Point 1 Longitude         [deg] (User)
%   h2         [1x1]       Point 2 Altitude          [km]  (Sat) [-180,+180]deg
%   lat2       [1x1]       Point 2 Latitude          [deg] (Sat) [-90,+90]deg
%   lon2       [1x1]       Point 2 Longitude         [deg] (Sat) [-180,+180]deg
%   data       [---]       Problem Data              [strc] (see NOTES)
%   f          [1x1]       Carrier Frequency         [Hz] (optional)
%--------------------------------------------------------------------------
% OUTPUTS:
%   TEC        [1x1]       Total Electron Content    [TECU]
%   drho       [1x1]       Range Error               [m]
%   tau        [1x1]       Ionospheric Time Delay    [ns]
%--------------------------------------------------------------------------
%
%NOTES:
% - The input "data" has been chosen to be a structure (for compactness of
%   the code) defined as:
%       data.a0      = 1st Az Coeff   [-]
%       data.a1      = 2nd Az Coeff   [-]
%       data.a2      = 3rd Az Coeff   [-]
%       data.mth     = Month          [month]
%       data.UT      = Universal Time [hours]
%       data.stModip = MODIP Table    [-] (modipneqg_wrapped.asc)
%       data.F2      = F2 Table       [-] (ccir21.asc)
%       data.Fm3     = Fm3 Table      [-] (ccir21.asc)
% - Giving in input the signal frequency the range error and time delay
%   will be computed.
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Initial Checks

%Check Input Data
[flag] = Checker(lat1, lon1, lat2, lon2, data);
if flag == 1
    return
end

%% Main Code

%Constants
RE = 6371.2;                %Earth Mean Radius [km]

%Ray-Perigee
[perigee] = RayPerigee(h1, lat1, lon1, h2, lat2, lon2);
%Radii of Objects [km]
r1 = h1 + RE;
r2 = h2 + RE;
%Intermediate Points [km]
s1 = sqrt(r1^2 - perigee.rp^2);
s2 = sqrt(r2^2 - perigee.rp^2);
%TEC Integration Strategy Determination
%Vertical TEC (P1 and P2 close)
if (abs(lat2 - lat1) < 1e-5) && (abs(lon2 - lon1) < 1e-5) %[2]
    %Vertical TEC (when P1 and P2 are too close)
    TEC = TECVert(h1, h2, perigee, lat1, lon1, data, 0.01);
else
    %Gauss Integration
    if (h1 < 1000) && (h2 > 2000)
        %Intermidiate Points
        sa = sqrt(54334589.44 - perigee.rp^2);
        sb = sqrt(70076989.44 - perigee.rp^2);
        %Integral Segmentation
        TEC1 = TECGauss(s1, sa, perigee, data, 0.001);
        TEC2 = TECGauss(sa, sb, perigee, data, 0.01);
        TEC3 = TECGauss(sb, s2, perigee, data, 0.01);
        %TEC
        TEC  = TEC1 + TEC2 + TEC3;
    else
        %TEC
        TEC  = TECGauss(s1, s2, perigee, data, 0.01);
    end
end

if nargin == 8
    %Time Delay [s]
    tau = 40.3/f^2 * TEC;
    %Speed of Light [m][s-1]
    c = 2.99792458e8;
    %Range Error [m]
    drho = tau * c;
end

end